Nonreciprocal forces enable cold-to-hot heat transfer between nanoparticles

We study the heat transfer between two nanoparticles held at different temperatures that interact through nonreciprocal forces, by combining molecular dynamics simulations with stochastic thermodynamics. Our simulations reveal that it is possible to construct nano refrigerators that generate a net heat transfer from a cold to a hot reservoir at the expense of power exerted by the nonreciprocal forces. Applying concepts from stochastic thermodynamics to a minimal underdamped Langevin model, we derive exact analytical expressions predictions for the fluctuations of work, heat, and efficiency, which reproduce thermodynamic quantities extracted from the molecular dynamics simulations. The theory only involves a single unknown parameter, namely an effective friction coefficient, which we estimate fitting the results of the molecular dynamics simulation to our theoretical predictions. Using this framework, we also establish design principles which identify the minimal amount of entropy production that is needed to achieve a certain amount of uncertainty in the power fluctuations of our nano refrigerator. Taken together, our results shed light on how the direction and fluctuations of heat flows in natural and artificial nano machines can be accurately quantified and controlled by using nonreciprocal forces.


A. Dragging experiment in MD simulations
To estimate the friction coefficient, a dragging MD simulation is performed assuming the Stokes' law for the friction force (F friction = −γv). In this regard, a constant force along the x axis of magnitude F = 0.3nN is applied to a spherical copper nanoparticle with radius r = 1.44nm immersed in an Ar fluid bath with fixed temperature 100K (first simulation) and 120K (second simulation) and periodic boundary conditions. As shown in Supplementary Fig.  S1, the velocity of the Cu nanoparticle reached a steady value (terminal velocity) after a transient time of ∼ 0.3ns and remained constant for the rest of the simulation. Using the Stokes' law, the friction coefficient is estimated from the ratio between the magnitude of the external force and the terminal velocity γ = 2.8 × 10 −12 kg/s obtained as the mean of the simulations done at the two temperature values. The simulation is repeated with an external force ten times larger (F = 3nN), yielding a very similar estimate of the friction coefficient γ = 3.1 × 10 −12 kg/s. Next, we describe how we extracted the effective friction coefficients from fits of the autocorrelation functions obtained from the MD simulations to analytical expressions derived below for our underdamped Langevin model. In thermal equilibrium, the normalised autocorrelation functions of the particle position and velocity are known. For κ C = κ H = 0, they read [1] and with ω 1 = κ m − γ 2 m 2 . Further, the position-velocity cross correlation is given by Here we calculate the ensemble-averaged heat flow rates ⟨Q j ⟩, with j ∈ {C, H}. To this end, we start from Sekimoto's definition [2] of the stochastic heat dissipated to the j bath in [t, t + dt], i.e. dQ j = (γẊ j − ξ j ) • dX j , and consider the steady-state average Inserting ⟨ξ jẊj ⟩ = γk B T j /m, which follows from the Langevin equation [3], one finds Thus, the heat flow between each nanoparticle and its bath is directly given by the variance of its velocity fluctuations. We note that this formula can also be recast into the form ⟨Q j ⟩ = kBγ m T eff j − T j , i.e., the heat flow is proportional to the difference between "effective temperature" T eff j = m⟨Ẋ 2 j ⟩/k B and actual temperature T j . The stationary joint probability density of the position and velocity of the two particles described by the linear Langevin equation [given in Eq (1) in the Main Text] is given by the four-variate normal distribution where is a column vector with T denoting matrix transposition, ⟨r⟩ its stationary average, and C the correlation matrix which is given by To calculate the variance of the velocity fluctuations, we develop further the formalism introduced in Refs. [4][5][6] to calculate correlation matrices. To this end, we first recast the Langevin equations into the matrix equationṡ with the state vector r = (X C , X H ,Ẋ C ,Ẋ H ) T and the noise vector Ξ = (0, 0, ξ C , ξ H ) T . The 4×4 matrix B contains the 2 × 2 zero and identity matrices, 0 2 and I 2 , the friction (diagonal) matrix Γ and the coupling matrix A. For the sake of generality, we will here consider the most general case of coupling here (four independent coupling entries of the coupling matrix), and below specialize our results to . Further, we define the diffusion From these ingredients, we can determine the correlation matrix C from the following expressions: where Q is an anti-symmetric 4 × 4 matrix that is uniquely determined by From this expression, and the fact that in stationary state with nonzero elements given by , (11e) This generalizes the results given in [5] to the case of two independent trap stiffness and coupling strengths. For the case considered here, Hence, we find substituting (12) in Eq. (5) which coincides with Eq. (3) in the Main Text. The expression (4) in the Main Text for the rate of heat dissipation to the hot bath can be found following analogous steps.
The last term denotes the rate of change of the Shannon entropy ⟨S sh ⟩ = k B ⟨− ln ρ 4 (X C , X H ,Ẋ C ,Ẋ H )⟩ of the joint probability density function. This term is constant in the steady state, thus, ⟨Ṡ sh ⟩ = 0. In the present case, using (3) and (4) from the Main Text, and ⟨Ṡ sh ⟩ = 0 and Eq. (14), we find which is also given in (5) in the Main Text. As expected, the total entropy production rate of the process is always greater or equal than zero. We further note that the rates of entropy production and heat dissipated have contributions from the inertial terms (i.e., it depends on m). Figure 3b in the Main Text shows the rate of entropy production from Eq. (15) and from the MD simulations, with the latter given by the heat flows extracted by the thermostat divided by the respective temperatures.
III. DETAILED BALANCE AT THE "PSEUDO EQUILIBRIUM" POINT As described in the Main Text, we find that at κ C T H = κ H T C the average values of the heat flows and the entropy production rate all vanish. To further investigate this "pseudo equilibrium" point, we check whether in this point the Langevin Eq. (8) fulfills detailed balance, meaning that all probability flows vanish. Detailed balance is a fundamental law defining of systems in thermal equilibrium. We employ a reasoning inspired by the arguments used in [8].
To this end, we consider the flow of the 4-point joint probability density function (pdf), ρ 4 (x, v, t), of x = (x C , x H ) T and v = (ẋ C ,ẋ H ) T , appearing in the corresponding multivariate Fokker-Planck equation. The Fokker-Planck equation reads with the probability currents J v , J x and the diffusion, friction and coupling matrices defined in (8). In general, the latter are constant in steady states, and zero in equilibrium. Now we use the identity ∂ z ρ = [∂ z ln(ρ)]ρ, and define the four-dimensional phase space velocity u to rewrite the Fokker-Planck equation (16) as The phase space velocity is connected to the probability current by J = In summary, this reasoning has revealed that the detailed balance is fulfilled if and only if (19), which coincides with the "equilibrium condition" found from the entropy production rate. We stress that this condition is irrespective of the friction coefficients, the trap stiffness and the mass of the nanoparticles, and has the same form when starting with the overdamped limit of the Langevin equations [9]. Note that the linear plots in (a,b) appears much more noisy, as few data points lie in the shown range, which is a consequence of the power-law distribution. The red lines in (c,d) depict power laws with exponent −2. The exponent of −2 is in good agreement with the COP of heat pump and of the refrigerator for all values of κH/κC considered.
In the Main Text, we consider the system as a nano refrigerator, i.e., we quantify the efficiency of cooling. However, one can also use this machine as a heat pump to heat up the hot bath. The efficiency of heat pumping can be measured by an analogous coefficient of performance defined as Like the COP of the refrigerator mode, the COP of the heat pump is bounded by the Carnot value ⟨ϵ p ⟩ ≤ TH TH−TC . Supplementary Fig. S3 displays both COPs for the refrigerator and heat pump, from simulations and from the theoretical predictions. Supplementary Fig. S4 shows the distributions of the fluctuating COPs, which are discussed in the Main Text. Here, we derive the distribution of the instantaneous (stochastic) power. We make use of the explicit expression for the stationary joint probability density for the position and velocities of the two nanoparticles ρ 4 (x C , x H , v C , v H ), which is a four-variate Gaussian distribution with zero mean and covariance matrix C (7) as given in this supplemental in Sec. II.
To access the distribution ofẆ j with j ∈ {C, H}, we first recall thaṫ with l ̸ = j which we assume throughout this section. From ρ 4 , we can derive P (Ẇ j ) by the following integral Here, is the joint stationary distribution of the position of the l−nanoparticle and the velocity of the j−nanoparticle, with j ̸ = l. It is given by a bivariate normal distribution with normalization constant N j , and coefficients which are functions of the covariances of the position and velocities of the nanoparticles. Therefore, the distribution of the power depends on the functions α j , β j , ζ j defined in (25), which in turn depend on ⟨X 2 l ⟩, ⟨X lẊj ⟩ and ⟨Ẋ 2 j ⟩. See Sec. II A of this supplemental for explicit analytical expressions of the covariances in terms of the physical parameters of the system.
Using the explicit expression of ρ 2 given by Eq. (24) and integrating out the delta distribution in Eq. (23), we obtain Changing variables y = √ α j x l and absorbing the Jacobian of the transformation in the new normalization constant Eq. (28) yields the distribution of the power exerted on the j ∈ {C, H} nanoparticle which coincides with Eq. (9) in the Main Text.
From the distributions (29), we can further deduce analytical expressions for any moment of the power. In particular, the variance of the power reads: for the cold particle (ẆC, blue), the hot particle (ẆH, red) and the total power (Ẇ =ẆC +ẆH, green). The values obtained from MD simulations are shown as symbols, the theoretical prediction are shown as lines.

VI. STABILITY CONDITIONS
Here we study the stability of the process described by the Langevin equation given in Eq. (1) in the Main Text, using the matrix equation given in (8). The stability conditions can be found from the roots of the eigenvalues of the generalised coupling matrix −B. The eigenvalues read The two largest eigenvalues are those with the respective plus sign. From them, we find the stability conditions −γ + γ 2 − 4mκ < 0, thus κ > 0, and −γ + γ 2 − 4m(κ + κ H + κ C ) < 0, thus κ H + κ C > −κ. To summarize, the system is stable if κ > 0, and κ H + κ C > −κ.
We note that these conditions are independent of m and coincide with the stability boundaries of the overdamped system (as expected). Further note that the eigenvalues become complex if γ 2 < 4mκ, or γ 2 < 4m(κ + κ C + κ H ). For the parameter choices considered in the MD simulations, all eigenvalues have imaginary parts, indicating oscillatory behaviour (for all considered κ H , κ C values). The oscillatory character of the stochastic dynamics (which stems from